%% Compare shock simulations in stationarised and unit-root models
% by Jaromir Benes

% In this m-file, we show that the three model versions have identical
% shock responses.

%% Clear workspace and load models

clear;
close all;
home;

%% Load Model Objects
%
% Load all three versions of the model created in `read_model`:
%
% * `m1` is the stationarised version,
% * `m2` is the unit-root version solved around a point on the
% balanced-growth path where productivity `A = 1`,
% * `m3` is the unit-root version solved around a point on the
% balanced-growth path where productivity `A = 2`.

load  read_model.mat m1 m2 m3;

%% Simulate a productivity shock
%
% Simulate the same shock in the three versions of the model, and compare
% the simulations. They are all identical.

d1 = sstatedb(m1,1:40);
d1.u(1) = 0.10;
s1 = simulate(m1,d1,1:40);
s1 = dbextend(d1,s1);

d2 = sstatedb(m2,1:40);
d2.u(1) = 0.10;
s2 = simulate(m2,d2,1:40);
s2 = dbextend(d2,s2);

d3 = sstatedb(m3,1:40);
d3.u(1) = 0.10;
s3 = simulate(m3,d3,1:40);
s3 = dbextend(d3,s3);

dbplot(s1,0:40,{'y','c','i','r','k','dA'},'tight',true);
ftitle(['Model m1 ',comment(m1)]);

dbplot(s2,0:40,{'Y/A','C/A','I/A','r','K/A','dA'},'tight',true);
ftitle(['Model m2 ',comment(m2)]);

dbplot(s3,0:40,{'Y/A','C/A','I/A','r','K/A','dA'},'tight',true);
ftitle(['Model m3 ',comment(m3)]);

disp([ s1.y, s2.Y/s2.A, s2.y, s3.Y/s3.A, s3.y ]);
disp([ s2.Y/s2.A-s1.y, s2.y-s1.y, s3.Y/s3.A-s1.y, s3.y-s1.y]);
maxabs([ s2.Y/s2.A-s1.y, s2.y-s1.y, s3.Y/s3.A-s1.y, s3.y-s1.y])

%% Unit-root responses
%
% With `m2` and `m3`, the responses of the unit-root variables (note the
% shock has permanent effect on all of them) can be obtained directly. This
% is convenient, isn't it?

dbplot(d2 & s2,0:40,{'Y','C','I','r','K','dA'},'tight',true);
legend('No-shock control','Shock simulation');
ftitle(['Model m2 ',comment(m2)]);

dbplot(d3 & s3,0:40,{'Y','C','I','r','K','dA'},'tight',true);
legend('No-shock control','Shock simulation');
ftitle(['Model m3 ',comment(m3)]);

%% Help on IRIS functions used in this m-file
%
% Use either `help` to display help in the command window, or `idoc`
% to display help in a HTML browser window.
%
%    help model/sstatedb
%    help model/simulate
%    help data/dbextend
%    help data/dbplot
%    help grfun/ftitle
